Optimal quantum control of Bose Einstein condensates in magnetic microtraps 
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Transport of Bose-Einstein condensates in magnetic microtraps, controllable by external param- 
eters such as wire currents or radio-frequency fields, is studied within the framework of optimal 
control theory (OCT). We derive from the Gross-Pitaevskii equation the optimality system for the 
OCT fields that allow to efficiently channel the condensate between given initial and desired states. 
For a variety of magnetic confinement potentials we study transport and wavefunction splitting of 
the condensate, and demonstrate that OCT allows to drastically outperfrom more simple schemes 
for the time variation of the microtrap control parameters. 



PACS numbers: 03.75.-b,39.20.+q,39.25.+k,02.60.Pn 



I. INTRODUCTION 

Trapping and coherent manipulation of cold neutral 
atoms in microtraps near surfaces of atomic chips is a 
promising approach towards full control of matter waves 
on small scales 0, H This field of atom optics is 
making rapid progress, driven both by the fundamen- 
tal interest in quantum systems and by the prospect of 
new devices based on quantum manipulations of neutral 
atoms. Lithographic and other surface-patterning pro- 
cesses nowadays allow to build complex atom chips which 
combine many traps, waveguides, and other elements, in 
order to realize controllable composite quantum systems 
[3| as needed, e.g., for the implementation of quantum 
information devices Such microstructured surfaces 
have been highly successful and form the basis of a grow- 
ing number of experiments [6(. 

The possibility to store, manipulate @, & S El El El , 
and measure a single quantum system with extremely 
high precision has initiated great stimulus in various 
fields of research, ranging from atom interferometry 
El El El El 51, over quantum sates [H, El [Tp 
and resonant condensate transport to microscopic 
magnetic- field imaging [2p| . In the vast majority of these 
schemes the wavefunction of the Bose-Einstein conden- 
sate, trapped in the vicinity of an atom chip, is manip- 
ulated through variation of the magnetic confinement 
potential. This is achieved by changing the currents 
through the gate wires mounted on the chip or mod- 
ifying the strength of additional radio-frequency fields 
d EH HI HI IHl ■ These external, time-dependent pa- 
rameters thus provide a versatile control for wavefunction 
manipulations, and make atom chips attractive candi- 
dates for quantum control applications. 

Consider the situation where one is aiming for an ef- 
ficient wavefunction transfer from a given initial to a fi- 
nal desired state, possibly by passage through a series 
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of other states, or for a conditional quantum gate where 
atoms have to interact with each other in a well-defined 
manner. Here the question arises: how should one mod- 
ify the control fields in order to achieve a most efficient 
transfer or coupling? This problem was first tackled by 
Hansel et al. [13| for a trapped-atom inteferometer setup 
where a dilute condensate should be split through varia- 
tion of the confinement potential from a single to a double 
well, such that it ends up in the groundstate of the final 
double well potential. These authors devised a scheme 
that optimizes adiabatic transfer by minimizing transi- 
tions to excited states. 

In this paper we study quantum control of Bose- 
Einstein condensates in magnetic microtraps within the 
framework of optimal control theory (OCT). Here, the 
objective of the control is quantified through a cost func- 
tion, which is then minimized subject to the condition 
that the time dynamics of the condensate is governed 
by the Gross-Pitaevskii equation 0, [25| . We will show 
that optimal control theory provides a versatile tool for 
determining efficient control strategies, and is applica- 
ble for realistic confinement potentials, one- and two- 
dimensional problems, and nonlinearities in the conden- 
sate dynamics. Optimal control theory is a mathematical 
device that allows for a general determination of efficient 
control strategies [H [2?J , and has found widespread ap- 
plications for, e.g., molecules [H [29[, atoms [H HU, or 
semiconductors [32[ . We believe that there are a number 
of reasons that render OCT ideal for quantum control 
of condensates in atom chips. First, it is only the cost 
function that determines the optimal control. There is 
no additional bias, such as, e.g., in the adiabatic scheme 
where scattering losses are minimized throughout the 
whole transfer process, and consequently OCT allows to 
explore a larger portion of the control space. In addition 
no knowledge of the stationary solutions of the Gross- 
Pitaevskii equation is required in OCT, contrary to the 
adiabatic scheme where ground and excited states must 
be determined for every control configuration. Since the 
optimal control corresponds to a minimum in the control 
space, the solutions are robust with respect to small flue- 
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tuations of the external parameters, which can never be 
avoided in real experiments. Finally, decoherence effects, 
that also play a role in atom chips [3, [H, [H, 
naturally incorporated into OCT calculations 



V(r, A = 1) at time T. Let 
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We have organized our paper as follows. In sec. ILT1 wc 
introduce to the realm of optimal quantum control, and 
derive the optimality system for condensate transport in 
atom chips. In sec. IIIII we present results for condensate 
splitting in simple and realistic confinement potentials. 
We demonstrate that our scheme is applicable for effec- 
tive one- and two-dimensional geometries, and for non- 
linearities in the condensate transport. Finally, in sec. HVI 
we summarize and draw some conclusions. 



II. THEORY 



We consider a coherent ensemble of Bosc-Einstcin con- 
densed atoms confined in a potential V(r, A(t)) produced 
by a magnetic microtrap. A(i) is a control parameter 
that describes the variation of the confining potential 
when changing the external parameters, such as currents 
through the microtrap wires or frequenc y a nd strength 
of additional radio- frequency fields [3, [H, [2l[ (for details 
see below). Through X(t) it is possible to manipulate 
the Bose-Einstein condensate, e.g., to split and reunite 
it by varying the potential from a single to a double well 
and vice versa. We assume that X(t) is a single- valued, 
real parameter, although different situations, e.g., micro- 
traps controlled by several parameters, could be treated 
equally well. In the following we shall assume for simplic- 
ity that A(t) can only take values between zero and one. 
The mean-field dynamics of the condensate is described 
by the Gross-Pitaevskii equation [24], |25| 



iip(r,t) 



V(r,\(t))+g\iP(r,t)\*)*(,(r,t) 



(1) 

with g a coupling constant related to the scattering 
length of the atoms. 1 Suppose that initially the system 
is in the groundstate ipo for the potential V(r,X = 0). 
Upon varying X(t) in the time interval t £ [0, T] from 
zero to one, the system will pass through a sequence of 
states and will end up in the final state ip{T). Within 
the field of quantum control one is usually seeking for an 
optimized time evolution of X(t) that allows to channel 
the system from the initial state ipo at time zero to a de- 
sired state ipd at final time T. In accordance to Ref. [HI], 
we assume to be the groundstate for the potential 



1 We set h = 1, and measure mass in units of the atom mass and 
length in units of micrometers. For 87 Rb atoms the time and 
energy scales are then given by 1.37 milliseconds and 5.58 nano 
Kelvins, respectively. 



J(V,A) = -(1 



\(M^(T))\ 



A(i) dt (2) 



be the cost function that rates how good the final 
state ip(T) matches the desired state tpd, with (u\v) — 
J dr u* (r)v(r) the usual inner product. 2 The first term 
on the right-hand side becomes zero for ip(T) = tpd and 
maximal if final and desired state do not overlap. 3 The 
second term on the right-hand side favours control fields 
X(t) with a smooth time variation and is needed to make 
the quantum control problem well posed. 7 is a weight- 
ing parameter that determines the importance of the 
two different control strategies of wavefunction match- 
ing and smooth control fields. We shall use small 7 val- 
ues throughout, such that the cost function J{tp, X) is 
dominated by the first term. The control problem under 
consideration thus becomes the minimization of the cost 
function J(ip, A) subject to the condition that ijj(t) fulfills 
the Gross-Pitaevskii equation ([TJ) . 

Within the field of optimal control theory (OCT) one 
uses Lagrange multipliers to turn this constrained min- 
imization problem into an unconstrained one. For this 
purpose we define the Lagrange function 

L(iP,p,X) = J(V>,A) 

+$e(p,ii>- (~V 2 + V x +g\^\ 2 )i\ , (3) 



with the abbreviation (u,v) = J Q dt J dru*v, and p(t) 
the Lagrange multiplier. We next utilize the fact that 
the Lagrange function has a saddle point at the mini- 
mum of J(ip,X), i.e. all three derivatives SL/6ip, SL/Sp, 
and 5L/SX must be zero. Performing usual functional 
derivatives in eq. ([3|) we obtain after some variational 
calculation the following optimality system 



ix/j 
if, = 
7 A = 



\v 2 + V x +g\i>\ 2 j <• Ha! 
-^ 2 + V x + 2g\^\ 2 Jp + gi; 2 p* (4b) 
-HeM^b). (4c) 



We assume that the wavfunction ip(r, i) is normalized to one. In 
comparison to, e.g., ref. |24| . where the wavefunction is normal- 
ized to the number of atoms N in the condensate, the nonlin- 
earity parameter g in eq. JTJ is therefore assumed to incorporate 
N. 

In contrast to the i||i/>(T) — tpd\\ 2 cost function used in ref. |27|1 , 
in expression J2} the final wavefunction has to match the desired 
one only up to a global phase e** 1 . This allows the substraction of 
constant values in the confinement potential V(r, A), as discussed 
in more detail in appendixfA] which proves particularly useful for 
magnetic confinement potentials with large energy offsets |21| . 
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which has to be solved together with the initial and ter- 
minal conditions 

V>(0) = ^0 (5a) 
w{T) = -{^{T))^ d (5b) 
A(0) = 0, A(T) = 1. (5c) 

The right-hand side of eq. (|5b|) follows from the func- 
tional derivative 5 J / Sip. Notice that while the state equa- 
tion (|4a[) with initial condition ip(0) — ipo evolves forward 
in time, the adjoint equation (|4b[) with terminal condi- 
tion (|5b[) is marching backwards. 4 The control equation 
(|4c]) determines the optimal control. 

In most cases of interest one is not able to directly 
guess X(t) such that eqs. (j4al -c) are simultaneously ful- 
filled, and one has to employ an iterative scheme. In 
this work we follow ref. [27| and formulate a numerical 
algorithm that solves the optimality system (|4"al -c) for 
given initial and desired configurations V'o an d V'rfi respec- 
tively. To solve this problem, we apply a gradient-type 
minimization algorithm, which, starting from an initial 
guess for X(t), determines a search direction for an im- 
proved control. We first solve eq. (|4a|) with initial condi- 
tion V'(O) = ipo forwards in time. Once the wavefunction 
ip(T) at time T is computed, the final condition for p(T) 
can be calculated from eq. (|5b[) and the adjoint equation 
of motion (|4bj) is solved backwards in time. The gradient 
of L with respect to A becomes 

^ = - 7 A-5Re(^|^| P >, (6) 

which gives the search direction for an improved control 
that minimizes J(tJj,X). In the following we employ for 
the minimum search either the usual gra dient method or 
the BFGS quasi-Newton method [27L l38j|. Details of our 
solution scheme for the Schrodinger-type equations are 
given in appendix [Bl 

We emphasize that the choice of our cost function ([2|) 
is by no means unique. For instance, one could add an 
additional 77 (ip, (— V 2 /2 + V\+ g\tp\ 2 ) tp) term, with an- 
other weighting parameter 77, to minimize the total en- 
ergy within the transfer process and to favour adiabatic 
processes. Another possibility would be to make 7 in 
eq. @ time dependent and to penalize control variations 
more strongly at the beginning and end of the transfer 
process, such that X(t) is turned on and off sufficiently 
smooth, which might be beneficial for experimental im- 
plementations. Alternatively, through a slight variant 
of the cost function the system can be forced to pass 



4 For the optimal control field \(t) the adjoint equation l|4b|l de- 
scribes the fluctuations of the system around ip(t). For this 
reason it is often referred to as the sensitivity equation. In- 
cidentally, eq. l|4b|l for the adjoint variable p closely resembles 
the time-dependent Bogoliubov— de Gennes equation [25j which 
usually serves as a starting point for the treatment of collective 
excitations. 
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FIG. 1: (Color online) Results of our simulation for the single- 
well potential (fTJl and xo = 5. (a) Cost function J(ip, A) for 
different transfer times T and for a linear time variation of A 
(solid line) and optimized variations (symbols). The dotted 
line shows results of a simulation where an additional x 4 term 
is added to the potential (see text). The inset reports the 
optimized control fields X(t) for transfer times of 3 (solid line), 
6 (dashed line), and 9 (dotted line), (b) Time evolution of 
potential and wavefunction \ip(x, t) \ for linear A and T=9. (c- 
e) Density plot of \ip(x, t)\ for linear A as a function of time 
and position, and for transfer times of (c) T=9, (d) T=6, and 
(c) T=3. On the bottom of each panel we show the desired 
wavefunction (dashed-dotted line) and \tp(x,T)\ (solid line). 
The solid lines in the density plots represent the equipotential 
lines of the confinement potential, (c'-e') Same as (c-e) but 
for optimized control fields. 

through a number of desired states [23] or to acquire a 
certain phase [17[ . One could also use a slight variant of 
our approach to obtain an optimization scheme for spa- 
tial geometries of waveguides and two-arm interferome- 
ters through which a condensate can propagate without 
creating excitations, as briefly outlined in appendix [U] 

III. RESULTS 

We first consider the more simple scenario of single 
atoms or dilute condensates within microtraps, and ne- 
glect the nonlinear terms in eqs. (T4"aj) and (|4b[) by setting 
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FIG. 2: (Color online) Time-integrated Wigner function 
w(x,p) for the single-well transfer processes shown in fig. [T] 
The solid lines show the equipotential lines for the Wigner 
functions of the initial and final wavefunctions ipo and ip{T), 
respectively. 



3 = 0. The influence of nonlinearities will be discussed at 
the end. In our optimal quantum control calculations we 
set 7 = 10~ 3 (smaller values of 7 turned out to have no 
noticeable influence on the results) and terminate the op- 
timization loop after several tens to hundreds iterations 
when the gradient ((6]) has become sufficiently small. 

We shall consider the scenario where a Bosc-Einstein 
condensate is split into two parts through smooth varia- 
tion of the magnetic confinement potential from a single 
to a double well [l(J, In sees. IIII Al and B we simu- 
late condensate transport through parabolic-like confine- 
ment potentials. These simplified case studies will allow 
us to grasp the essentials of our optimal control calcula- 
tions. Transport through realistic magnetic confinement 
potentials is discussed in sec. IIII CI and D. Finally, in 
sec. IIII El we investigate the influence of nonlinearities in 
the Gross-Pitaevskii equation. 



A. Single well 

After successful splitting of the condensate the atoms 
in the two respective wells can be further transported by 
shifting the location of the minima. In our first example 
we will study such transport inside a single well. We 
will make the assumption that the confinement along y 
and z is much stronger than along x, such that only the 
dynamics in x direction has to be considered. We assume 
a potential of the form 



which has its minimum at Xxq. By varying within the 
time interval [0, T] the control parameter from zero to 
one, the potential minimum becomes shifted from zero 
to xq. Our objective now is to seek for a time varia- 
tion X(t) that brings the system from the groundstate ipo 
of the harmonic oscillator centered at x — to the de- 
sired groundstate tpd of the displaced harmonic oscillator 
centered at xq- Although the above model allows un- 
der quite general conditions for an analytic solution, see, 
e.g., Ref. [23], the following analysis will turn out to be 
helpful when discussing the more complicated situation 
of condensate splitting. 

Let us first consider a linear variation X(t) = t/T. 
Figs. QJb-e) report results of our simulations for three 
selected tranfer times T. Panel (b) shows for T = 9 the 
modulus of the wavefunction together with the confine- 
ment potential at different times, and panel (c) shows a 
density plot for |-0(a;, i) | of the same transfer process. At 
the bottom of panel (c) we also plot the final wavefunc- 
tion (solid line) which somewhat differs from the desired 
one (dashed-dotted line). Similar behavior is observed 
for (d) T = 6 and (e) T = 3. Finally, in panel (a) we 
report for the linear A- variation the cost function © as 
a function of T (solid line) , which is high for small values 
of T and shows an oscillatory behavior with decreasing 
amplitude for longer transfer times T . The decreasing 
amplitude is due to the fact that with increasing T the 
time variation of potential ([7]) becomes slower, and the 
system can follow almost adiabatically. 

The oscillations in the cost function are due to the os- 
cillations of the wavefunction inside the single-well po- 
tential. To understand their origin, consider the ex- 
treme case where the position of the potential minimum 
is abruptly moved to x$ at time zero, and the system 
is brought into a highly excited state where the ground- 
state wavefunction of the harmonic oscillator is displaced 
by xq with respect to the new minimum of V(x, A = 1). 
Such displaced groundstates of the harmonic oscillator 
are known as coherent states [40j and have a dynamics 
reminiscent of classical oscillators. As time goes on, the 
system will start to oscillate with amplitude xq around 
its new equilibrium position. Also the A(i) variation with 
finite speed can be described in terms of such coherent 
states, as evidenced by the fact that in figs. DJb-e) only 
the position but not the shape of the wavepacket changes 
with time. 5 To inquire more into this evolution we shall 
analyze the Wigner function [411 ] 



w(x,p;t) 



e~ lps i>(x + -,t)ip*(x - -,t)ds (8) 



1 2 

V{x, A) = - (x ~ \x ) 



(7) 



5 This situation corresponds to the classical analogon of a particle 
attached to a spring which is initially fixed at the origin, and the 
point where it is fixed is moved to a different position xq at later 
time. Only for certain transfer times T the particle will end up 
in rest position and with no force acting upon it. These times 
correspond to the minima of the cost function in fig. [Ha). 
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FIG. 3: (Color online) Same as fig. [T] but for the double- well 
potential {SJ. The inset in panel (a) shows the optimized X(t) 
for T = 6 (solid line) and T = 9 (dashed line), (b) Time 
evolution of |^>(a;,i)| and V(x,X(t)) for T = 9. (c,d) Density 
plot of \ip(x,t) \ for (c) T = 9 and (d) T = 6. (b'-d') Same as 
(b-d) but for optimized control fields. 



for the wavefunction, which is a mixed position- 
momentum distribution that has many, albeit not all, 
properties of a classical distribution function. The solid 
lines in fig. [2] show contour lines of the Wigner functions 
for the initial and final wavefunctions i/jq(x) and ip(x, T), 
respectively. These coherent states are minimum uncer- 
tainty states with Ax Ap — | . The density plots in the 
different panels of the figure show the time integrated 
Wigner function w(x,p] t) dt which provides informa- 
tion about the trajectory in phase space. At short trans- 
fer times, fig. EJc), the system ends up in a state that is 
located close to xq but with a high momentum. Thus, 
when the control parameter is kept fixed to A = 1 at 
times beyond T, the system will continue to oscillate 
around its new equilibrium position. This final state dif- 
fers substantially from the desired groundstate ipd of the 
displaced oscillator, and consequently has a rather high 
cost J [see also fig.QJa)]. With increasing T [figs, [^a, b)] 
the momentum of the coherent state decreases and thus 
J becomes smaller. 

We next turn to our optimal quantum control calcula- 
tions. Here we start from the linear X(t) = t/T function 
as a guess for the control, and successively improve X(t) 
according to the scheme described in the previous sec- 
tion. The symbols in fig. [lja) report the cost function 
for the optimized control fields: throughout J can be 



drastically improved with respect to the linear A varia- 
tion. In particular for transfer times beyond say T = 3 
the final wavefunctions perfectly match the desired one. 
Figs. QJb'-e') report the wavefunction evolution for the 
optimized control fields depicted in the inset of panel (a). 
For T — 6 and T = 9 the fields deviate only little from the 
linear dependence, and minor to moderate corrections of 
X(t) suffice to finally bring the system at Xq to rest. This 
is also apparent from the Wigner functions shown in fig. [2] 
where the final state (solid contour line) is centered at xq 
and has zero momentum. 6 For T — 3 the control strat- 
egy shown in panel (e') becomes noticeably modified: at 
early times the center of the parabolic confinement po- 
tential is quickly shifted and the system is put into a 
highly nonequilibrium state where it starts to oscillate 
from left to right. Once it decelerates and reaches the 
right turning point, the position of the potential mini- 
mum is further shifted and the system becomes frozen in 
the groundstate of the shifted harmonic oscillator. Note 
that only the quasi-Newton BFGS method is capable of 
coming up with such control, whereas the more simple 
gradient scheme cannot and is trapped in a sub-optimal 
extremum. 

We emphasize that all optimal transport processes dis- 
cussed here rely on non-equilibrium coherent states, and 
the resulting transfer processes strongly differ from adi- 
abatic schemes. Finally, in panel (a) the dotted line re- 
ports that similar behavior is also found when an addi- 
tional nonlinear potential term 7j(x — Xxq) /4 is added, 
with rj = 0.2 in the figure. Also in this case optimal con- 
trol theory gives control fields (not shown) that allow, in 
contrast to the linear A variation, perfect transport. 

B. Double well 

We next turn to the more complicated situation of 
wavefunction splitting. As a preliminary case study we 
consider the confinement potential 

fi(M-f) 2 for 1*1 >f 
V(x,X) = l a (9) 
I I a;2 ) otherwise, 

which changes from a single- well potential for A = to a 
double- well potential with interwell distance d for A = 1 . 



6 It is worth noting that our cost function (J2J) is only governed 
by the final wavefunction ip(T), and consequently no guidance 
of the intermediate wavefunction trajectory is present. Thus, 
if the linear control fields already work successfully, such as for 
T = 6 in fig. [2d), the optimized A(t) and the corresponding 
transfer process are practically not altered, whereas somewhat 
stronger deviations can be observed for T = 9 in fig- [He). We 
also emphasize again that the parabolic confinement potential 
171 is special in the sense that it can only change the position 
but not the shape of the initial wavepacket, and there thus exists 
a huge variety of different successful control stratgies. 
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FIG. 4: (Color online) Absolute value of Wigner function 
\w(x,p; T)\ at the end of the transfer process [figs. EJb-d)] for 
linear A variation and for transfer times of (a) T = 9 and (b) 
T = 6. (a',b') Same as (a,b) but for optimized control fields. 
The interference pattern around the origin is a measure of the 
coherence properties of the split condensate. 



Potential © is constructed such that it is continous and 
smooth. Figs. EIJb,c) show results for a wavefunction 
splitting for T = 9 and for a linear \(t) dependence. 
The wavefunction becomes split in the first stage of the 
time evolution, and is transported in the respective min- 
ima in the second stage of the transport process. Con- 
trary to the single- well transport, in this second stage 
also excited vibrational states of the harmonic oscilla- 
tor that were populated during the initial splitting pro- 
cess are involved, as apparent from the varying shape 
of |^>(a;,£)| in the density plot of panel (c). Even more 
striking, the wavefunction shown in panel (d) for the fast 
transfer process with T — 6 is split only incompletely, 
and part of the population remains localized in-between 
the two wells. Correspondingly, the overlap with the de- 
sired groundstate ipd of the V(x, A = 1) potential is rather 
poor and the cost function shown in fig. 03a) (solid line) 
is high for small values of T. With increasing T the 
cost function again exhibits an oscillatory behavior with 
decreasing amplitude, indicating the onset of adiabatic 
transport. However, in contrast to the single- well case 
J keeps a finite value at its minima, which is due to the 
population of excited vibronic states during splitting and 
the resulting lack of complete overlap with t/^. 

The symbols in fig. (Ha) report the cost function for 
the optimized process of wavefunction splitting. At short 
transfer times, say below T — 5, the optimized con- 
trol strategies perform significantly better than the linear 
ones, but the overlap with the desired state is not per- 
fect. We emphasize that these results do not exclude 
the possibility of more efficient transfer in regions of the 
control space that were not explored by our minimiza- 
tion scheme. For transfer times beyond T = 6 the cost 
function drops below a value of 1CT 3 indicating the on- 
set of efficient wavefunction splitting. Figs. [3]Jb',c') show 
the wavefunction evolution for T = 9, which is not drasti- 
cally altered in comparison to the evolution for the linear 
scheme. A slight modification of the control function X(t) 



suffices to channel the system to the desired state at time 
T. Fig. 2] shows the Wigner functions w(x,p;T) at the 
end of the transfer processes. For the optimized control 
shown in panel (a') it consists of two coherent-state fea- 
tures at the positions of the two minima of the double 
well potential (see cross symbols), indicating that i/)(T) 
matches the respective single-well groundstatcs, and an 
interference pattern at position zero due to the superpo- 
sition nature of the wavefunction ip(T) [41(. In contrast, 
the Wigner function for the linear time evolution shown 
in panel (a) exhibits an asymmetric shape at the posi- 
tions of the potential minima, that can be traced back to 
the superposition of ground and excited vibronic states 
within the respective minima. For the short transfer time 
of T = 6 the optimized X(t) shown in the inset of fig.[3^a) 
substantially differs from a linear behavior. As apparent 
from the corresponding wavefunction evolution shown in 
fig- Eld'), at early times the potential is quickly trans- 
formed from a single to a double well, and the system 
is thereby brought into a highly excited state. Similar 
to the fast single-well transport described above, such 
states can be manipulated and transported on shorter 
timescales. Indeed, the final stage of the transfer process 
is reminiscent of the final stage of wavefunction transport 
shown in fig. Llje'). 

C. Magnetic confinement of Hansel et al. 

In Hansel et al. [l3| the authors studied wavefunction 
splitting for a realistic magnetic microtrap. They de- 
vised a control scheme that favours adiabatic transport 
by minimizing, throughout the whole transfer process, 
excitations to excited states, and demonstrated that such 
approach can perform significantly better in comparison 
to more simplified control strategies. In this section we 
re-examine their scheme within the framework of opti- 
mal control theory. We use the same model parameters 
for the magnetic microtrap 7 where confinement along x 
is provided by three parallel wires oriented along the y 
direction, with an inter-wire distance of 20 /im. The cur- 
rent J ext through the central wire is opposite to the cur- 
rents I c through the outer wires. Introducing a current 
modulation by means of the control parameter A via 

/ cx t = 140 + A x 2.91 mA 
I c = 0.25 + A x 4.4 mA (10) 

produces a magnetic confinement along x that changes 
from a single well at A = to a double well at A = 1, 
as shown in fig. [5fb) . For the linear variation of \(t) 
wavefunction splitting is shown in panels (c) and (d) for 
transfer times of 15 and 8 ms, respectively. In both cases 



7 The strength of the field component Bq,„ should be 30 G, rather 
than the 20 G given in eq. (2) of ref. |l3j . in order to match the 
distance of 35 (im between trap and surface. 
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FIG. 5: (Color online) Results for the magnetic microtrap of 
Hansel et al. [l3| and for 87 Rb atoms confined in the mp — 2 
state, (a) Cost function for linear t/T (solid line), square- 
root y/tjT (dashed line), and optimized (symbols) variation 
of The optimized J (if), A) is magnified in the left in- 

set. The right inset reports the optimal A(t) for T = 8 ms 
(solid line) and T — 15 ms (dashed line), (b) Contour plot 
of magnetic confinement potential as a function of A. (c,d) 
Wavefunction evolution for linear variation of A and transfer 
times of (c) T = 15 ms and (d) T = 8 ms. (c',d') Same as 
(c,d) but for optimized control. 



the splitting is too fast to allow the system to become 
localized in the two mimima of the double well, and a 
significant portion of the population remains in-between 
the two wells. This is also apparent from the cost func- 
tion shown in fig. [^a) (solid line) that reports large J 
values over a wide range of transfer times, thus indicat- 
ing an only incomplete splitting. The relation of our cost 
function to the excitation probability p used in ref. [l3| 
is simply given by J(ip, A) ~ hp, assuming as usual only 
minor contributions from the second term in eq. @. 

The symbols in fig.EJa) show that optimal control the- 
ory again allows to strongly improve the cost function. 
In the inset we report that for transfer times beyond say 
6 ms the cost function J becomes significantly lower than 
the control penalization 7 = 10~ 3 (dotted line), and the 
final wavefunction ip(t) matches almost perfectly the de- 
sired groundstate wavefunction of the final double-well 
potential. A comparison of the optimal control X(t) de- 
picted in the second inset of panel (a) with the optimized 
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FIG. 6: (Color online) Results for the magnetic microtrap 
of Lesanovsky et al. [5l|] and for 87 Rb atoms confined in 
the rriF = 2 state, (a) Magnetic confinement potential for 
different rf field strengths, (b) Linear (dashed line) and opti- 
mized control parameter A(i) as obtained from the solutions of 
the one- dimensional (solid line) and two-dimensional (dotted 
line, indistinguishable from solid line) Schrodinger equation, 
(c) Density plot of wavefunction evolution J \ip(x,y; t)\ dy for 
linear A variation, (d) Same as (c) but for optimized \{t). 
(e) Density plot of wavefunction evolution J \ij)(x, y, i)\ dx for 
optimized A(i). 



control of Hansel et al. , see inset of fig. 6 of ref. [l3[ , shows 
that both control strategies are of equal simplicity. We 
note that the optimal control fields of our approach per- 
form better for very short transfer times, whereas for 
longer transfer times further analysis would be needed to 
pinpoint the advantages and disadvantages of the respec- 
tive schemes. 



D. Magnetic confinement of Lesanovsky et al. 

In our fourth case study we consider the radio- 
frequency double-well confinement proposed by 
Lesanovsky et al. [2l| which is produced by a 
surface-mounted dc four-wire structure on an atom 
chip. Such traps provide tight confinement even at 
large surface distances, allow for smooth potential 
transitions by variation of external parameters, such 
as rf field strengths, and are relatively robust against 
experimental fluctuations. In our calculations we use the 
same parameters as given in ref. [2l| (see also eq. (10) 
therein) , and vary the rf field strength by means of the 
control parameter A according to 



B lf = 0.5 + A x 0.3 G. 



(11) 



Fig. Ela) shows the confinement along x for three dif- 
ferent rf field strengths corresponding to A = (solid 
line), A = A (dashed line), and A = 1 (dashed-dotted 
line). Contrary to the double- well potential discussed 
in the previous section, the magnetic confinement of 
Lesanovsky and coworkers exhibits a substantial exten- 
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sion in y direction, which calls for a solution of the two- 
dimensional Schrodinger equation. In our work this is 
accomplished by using the split operator technique, as 
discussed in more detail in appendix [BJ Figs. [Sfb-e) 
show results of our optimal quantum control calculations 
for a quite short transfer time of T = 2 ms. In panel (b) 
we report the optimized A(t) functions as obtained from 
the solutions of the one-dimensional (solid line) and two- 
dimensional (dotted line) Schrodinger equation. Both 
A(i) arc almost identical. Indeed, from fig. [He), which 
shows the wavefunction and confinement potential along 
y, it is apparent that there is an only minor influence 
of B T { on the confinement along y, and consequently the 
wavefunction factorizes. The lower parts of panels (c-e) 
report the final (solid line) and desired (dashed-dotted 
line) wavefunctions. They differ in case of a linear varia- 
tion of X(t) [see panel (c)] and coincide for the optimized 
control [see panels (d,e)]. Thus, optimal quantum con- 
trol allows to drastically outperform more simple con- 
trol schemes. We have presented the results of fig. [B] 
primarily to demonstrate our ability to also cope with 
two-dimensional problems. We believe that this will be 
important for the future analysis of more complicated 
potentials, such as ring-shaped interferometers [2l| . In 
our concluding remarks we will further elaborate on this 
point. 



E. Solution of nonlinear Gross— Pitaevskii equation 

Let us finally address the influence of the non-linear 
term in the Gross-Pitaevskii equation {T]) on our optimal 
quantum control results. In doing so we shall re-examine 
the results of sec. IIIIBI for the simple double dot po- 
tential ([9]), though similar results are also found for the 
more realistic potentials studied in sees. IIII Cl and UlI Dl 
Consider the one-dimensional Gross-Pitaevskii equation 

fid 2 \ 
iiP(x,t) = l--— + V(x,X(t)) + K\^(x,t)\ 2 \ i/>(x,t), 

(12) 

where all details of the condensate density and the 
transversal confinement potential have been lumped into 
the single nonlinearity parameter k. Figure [7]^ a) shows 
the groundstate wavefunctions ipo (%) (see appendix IB 21 
for computational details) and the effective potentials 
V e ff(x, A) = V(x, X) + n\ip (x)\ 2 for a few selected n values 
and for A = 0. Due to the repulsion of atoms in the con- 
densate the wavefunction broadens and penetrates into 
the barrier. 

We first consider condensate splitting through linear 
variation of A. Fig. [7jb) shows a density plot of the 
corresponding cost function J{ip, A) for different trans- 
fer times T and nonlinearity parameters k. Note that for 
k = the cost function J corresponds to the solid line 
shown in fig. [3Ja). From the figure we observe that for 
small transfer times, say below T = 5, the condensate 
becomes splitted only very inefficiently and there is no 
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FIG. 7: (Color online) Results of our simulations performed 
with the nonlinear Gross-Pitaevskii equation for the double- 
well potential ([9}. (a) Groundstate wavefunction f/>o and ef- 
fective potential Vo + ft|V'o| 2 for A = (single well) and differ- 
ent values of the nonlinearity parameter ft. (b) Density plot 
of the cost function J(tp, A) as a function of transfer time T 
and nonlinearity parameter k, and for a linear A variation. 
The contour at n = corresponds to the solid line shown in 
fig- EI a)- ( c ) Cost function as a function of n and for transfer 
time T = 8. The solid line and symbols correspond to the lin- 
ear and optimized variation of A, respectively. (d,d') Density 
plot of \ip(x, i) | for ft = 20 and T — 8, and for a (d) linear and 
(d') optimized A-variation. 



substantial overlap of the final wavefunction with the de- 
sired one. With increasing T the transfer process works 
more efficiently. Generally speaking, for comparable val- 
ues of J larger nonlinearities k translate to longer transfer 
times. This is also apparent from fig. [7]Jc) that reports 
J as a function of k for a fixed transfer time T = 8. 
The symbols in the figure show results of our optimal 
control calculations, based on the solutions of eqs. 
c) , which demonstrate perfect condensate splitting within 
the whole k regime under consideration. The correspond- 
ing time evolutions of the control parameters X(t) (not 
shown) are similar to those shown for k — in the inset 



9 



of fig. 02a). Thus, optimal control theory allows to devise 
efficient control strategies even in the presence of moder- 
ate condensate nonlinearities. Although the nonlinearity 
parameter influences the detailed time evolution of A(i), 
it has no drastic impact on the essentials and qualitative 
features of our findings. Similar conclusions apply to the 
results for other magnetic confinement potentials. 

In this paper we have only considered the mean- 
field Gross-Pitaevskii dynamics and have neglected 
Bogoliubov-type quasiparticle excitations out of the con- 
densate. This approximation is justified when the de- 
pletion of the initial Gross-Pitaevskii ground state due 
to quasiparticle excitations is sufficiently small (2f|. For 
given trap parameters, condensate density, and tempera- 
ture, this can in principle be determined from the noncon- 
densate normal and anomalous density matrices n and to, 
respectively, to be computed from the generalized Gross- 
Pitaeveskii equations [25|, |42|, 0, |44[. If h and rh are 
initially negligible, they will remain negligible through- 
out the transfer process. Otherwise the dynamics of the 
Bogoliubov quasiparticles should be explicitly accounted 
for [43l . l44j , which, within the framework of optimal con- 
trol theory, could be done by introducing additional La- 
grange parameters for n and to. Although such analy- 
sis is beyond the scope of the present paper, we expect 
from related studies |32l . [3?J that for this extended sys- 
tem OCT will not only allow to control the condensate 
wavefunction but also its quasiparticle excitations. 



IV. SUMMARY AND DISCUSSION 

In conclusion, we have studied within the framework of 
optimal control theory and the Gross-Pitaevskii equation 
quantum control of Bose-Einstein condensates in mag- 
netic microtraps, which can be controlled by external pa- 
rameters such as wire currents or radio- frequency fields. 
For a variety of magnetic confinement potentials trans- 
port and wavefunction splitting of the condensate has 
been analyzed, and we have demonstrated that OCT can 
drastically outperform more simple control strategies. 

In contrast to adiabatic transfer schemes, where the 
control fields in the time-dependent Hamiltonian H{t) 
must be changed sufficiently smooth, such that transi- 
tions to excited states are suppressed throughout, 8 OCT 
allows to access excited states during the control. This 
opens the possibility to explore a larger portion of the 
control space, and enables high-fidelity quantum control 
even at short time scales. Furthermore, neither the wave- 
functions nor energies of ground and excited states are 
needed in OCT calculations, which appears to be par- 
ticularly advantageous for microtraps controlled by sev- 



eral external parameters, where the solution of the time- 
independent Schrodinger or Gross-Pitaevskii equation for 
every configuration of the magnetic confinement poten- 
tial would be a computationally heavy task. OCT calcu- 
lations for Bose Einstein condensates in magnetic micro- 
traps are expected to be a useful and versatile tool for 
high-fidelity quantum control in a variety of applications. 
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APPENDIX A: GLOBAL PHASE 



Let 



V{r,\) = V {\)+V{r,\) 



(Al) 



where Vb(A) is the minimum value of the potential 
V{r, A) and V{r, A) a potential with minimum zero. 
Then, 



1>{t) = exp (-#(*)) i>(t) 



(A2) 



is a modified wavefunction with the global phase <f>(t) 
defined through 



*(*) = f V a 
Jo 



(X(s))ds. 



(A3) 



Inserting wavefunction (|A2[) into the Gross-Pitaevskii 
equation |T]) gives for the time derivative 

i<tP = e -^(*) (ijj + <j)(t) <$\ = ie-^ w 4 + V (X(t)) ip . 

(A4) 

Here, the last term on the right-hand side is canceled by 
a corresponding term on the right-hand side of eq. fT}, 
and we obtain the modified Gross-Pitaevskii equation 

4(r, t) = (-\^ 2 + V (r, A(t)) + g |$(r, t) |j $(r, t) . 

(A5) 

It is identical to eq. |T]) with the only exception that the 
constant contribution of the confinement potential V(t) 
is substracted. The solution ip(t) equals tp(t) up to the 
global phase <fr(i). 



APPENDIX B: DETAILS OF OUR NUMERICAL 
SCHEME 



1. Schrodinger equation 

More specifically, the inequality \(f\dH(t)/dt\0)\ < (E f - E ) 2 

must be fulfilled throughout the transfer process, with and / In this appendix we present details of our numerical 

denoting ground and excited state, respectively. solution schemes for tp(t) andp(i). Let US first neglect the 
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nonlinear terms in eqs. (|4a[) and (|4bp. by setting g = 0, 
and consider one spatial dimension. We discretize the 
time and space domain into a finite number N t and N x 
of subintervals of sizes St and Sx, respectively. A discrete 
state variable at time t m and position Xi is denoted by 
tpf 1 . The second spatial derivative t/j" is approximated by 
the finite difference expression (ipi+i ~ 2^ + ipi—i)/ (Sx) 2 
together with periodic boundary conditions. For the time 
integration of ip(t) and p(t) we use the Crank-Nicholson 
scheme 



1> 



rn+1 . 



= (1 + ^^i^ m+1 )- 1 (1 - i%R m )^ m . (Bl) 



The inversion of the matrix on the right-hand side is 
computationally simple owing to its tridiagonal shape, 
which results from the above finite difference scheme. 
The Crank-Nicholson scheme has the advantage that the 
time evolution is unitary and the norm of the wavefunc- 
tion is thus conserved. In our calculations we typically 
use values of N t = 500 and N x = 500. 

For two space dimensions we again discretize the do- 
main into equidistant subintervals and approximate the 
spatial derivatives by a corresponding five-point formula 
for the Laplacian. However, the resulting Hamiltonian 
matrix H is no longer tridiagonal and the inversion in 
(IB1[) becomes computationally more costly. We thus em- 
ploy the split-operator scheme [45j |. Let H = T + V, 
where T is the kinetic term resulting from the discretiza- 
tion of the Laplacian and V the magnetic confinement 
potential. The wavefunction at later time is then com- 
puted according to 



m+1 



iStT , 



(B2) 



which again is a norm-conserving scheme. Here, the ac- 
tion of the matrix e~ lStT on the wavefunction can be eas- 
ily computed by means of fast Fourier transform (FFT) 
and its inverse. Space discretizations of typical dimension 
512 x 128 can be easily handled within such approach. 

Finally, the second time derivative of A needed in 
eq. ([S]) is approximated by the finite difference expres- 
sion (A m+1 - 2A m + A m - 1 )/((5t) 2 . In our calculations we 
set 7 = 10~ 3 . 



2. Gross— Pitaevskii equation 

The split-operator technique can be also applied to the 
nonlinear Gross-Pitaevskii equation ([T]) . We first replace 
in eq. (|B2|) the potentials V by the effective potentials 
V c fi = V + g | %/} | 2 . An apparent difficulty of the nonlinear 
term is the fact that the first exponetial on the right- 
hand side of eq. (|B2|) invokes the wavefunction ip m+1 , 
through the effective potential V^2 + , which is not at 
hand at this stage of computation. One can, however, 



easily check that the e 1 ? 



term only adds a phase to 



modulus can be computed from 

.«<•-! I - e -i8tT e -i%VZ 



(B3) 



Once |?/> m+1 | is known we can determine the effective 
potential VJ% + , and finally compute the wavefunction 
at later time through 



m+l _ e ~^K 



1 'e- l5tT e- l ^«r 



(B4) 



This scheme again conserves the norm. As for the adjoint 
equation (|4b[) . we use a slight variant of the split-operator 
technique for the linear Schrodinger equation, where the 
real and imaginary parts of the equation are separated 
to cope with the g ijj 2 p* term. 

Finally, for imaginary time steps —iSt we are able 
to compute the groundstate wavefunction of the Gross- 
Pitaevskii equation. Here, we start from the groundstate 
of the linear Schrodinger equation and evolve the system 
through eqs. (|B3IB4|) in imaginary time, thus projecting 
out the groundstate wavefunction. After each iteration 
the wavefunction is normalized and the computation ter- 
minates when the wavefunction does no longer change 
significantly. 



APPENDIX C: OPTIMIZATION OF SPATIAL 
GEOMETRIES 

In this appendix we briefly discuss how a slight vari- 
ant of the OCT scheme presented in sec. |H] would allow 
for an optimization of spatial geometries, such as waveg- 
uides or two-arm interferometers. The situation we have 
in mind is a scattering- type experiment, where the con- 
densate, initially in state i/'O) enters through a waveguide 
into the scattering region where it becomes split. Let us 
consider for simplicity a two-dimensional geometry and 
a condensate propagation along x. The objective of our 
optimization thus becomes the choice of the confinement 
potential V(x,y) in the scattering region through which 
one can propagate the condensate without creating exci- 
tations. 

Let ipd denote the desired outgoing state of the scat- 
tering and V(y,X(x)) the confinement potential param- 
eterized through the space-dependent control parameter 
A (a;), with x s [0, L\, Instead of eq. @ we introduce the 
cost function 



JW>,A) = ~(1 



,(C1) 

where the last term favours a smooth spatial variation 
of the confinement potential V . Performing functional 
derivatives of the Lagrange function we obtain again 
eqs. (0Ji,b), whereas (0c]) has to be replaced by the space- 
dependent version 



the different components ip^ 1 + ■ Thus, the wavefunction 



dx 2 y\d\(x) 



(C2) 
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with boundary conditions A(0) = and A(L) = 1. Here solution of the resulting optimality system can be per- 
the expression on the right-hand side involves an into- formed along the same lines as for its time-dependent 
gration over the transversal coordinate y and time. The counterpart. 
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